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Existing theoretical models of evolution focus on the relative fitness advantages of different mu- 
tants in a population while the dynamic behavior of the population size is mostly left unconsidered. 
We here present a generic stochastic model which combines the growth dynamics of the popula- 
tion and its internal evolution. Our model thereby accounts for the fact that both evolutionary 
and growth dynamics are based on individual reproduction events and hence are highly coupled 
and stochastic in nature. We exemplify our approach by studying the dilemma of cooperation in 
growing populations and show that genuinely stochastic events can ease the dilemma by leading to 
a transient but robust increase in cooperation. 
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Commonly, Darwinian evolution in terms of reproduc- 
tion, selection, and variation is described in frameworks 
of population genetics and evolutionary game theory [T]- 
3 . These approaches model the internal evolutionary 
dynamics of a species' different strategies (or traits) in a 
relative perspective. Namely, they compare fitness terms 
and focus on the relative advantage and abundance of dif- 
ferent traits. In such a setup, the time evolution of the 
relative abundance, x, of a certain strategy is frequently 
described by a replicator equation, 



d t x = (f-(f))x. 



(1) 



A trait's relative abundance will increase if its fitness / 
exceeds the average fitness (/) in the population. 

While in these evolutionary approaches, the dynamics 
of the population size, N, is mostly left unconsidered or 
assumed to be fixed [3j, in population ecology the dy- 
namical behavior of a species' population size is studied. 
Models of population dynamics [H[S] usually describe the 
time development of the total number of individuals, N, 
by equations of the form 



d t N = F{N,t). 



(2) 



T (N, t) is in general a non-linear function which includes 
the influence of the environment on the population, such 
as the impact of restricted resources or the presence of 
other species. By explicitly depending on time a changing 
environment such as, for example, the seasonal variation 
of resources can be taken into account. 

The internal evolution of different traits and the dy- 
namics of a species' population size are, however, not 
independent [B|. Actually, species typical coevolve with 
other species in a changing environment and a separate 
description of both evolutionary and population dynam- 
ics is in general not appropriate. Not only population dy- 
namics affects the internal evolution (as considered, for 
example, by models of density-dependent selection [7]), 
but also vice versa. Illustrative examples of the cou- 
pling are biofilms which permanently grow and shrink. 



In these microbial structures diverse strains live, inter- 
act, and outcompete each other while simultaneously af- 
fecting the population size [5]. So far, specific examples 
of this coupling have been considered by deterministic 
approaches only, e.g. jHHU]. However, classical and re- 
cent work has emphasized the importance of fluctuations 
for internal evolution which are only accounted for by 
stochastic, individual-based models, e.g. ;2j 3l [T2| ITS"]. 

In this Letter, we introduce a class of stochastic models 
which consider the interplay between population growth 
and its internal dynamics. Both processes are based 
on reproduction events. A proper combined description 
should therefore be solely based on isolated birth and 
death events. Such an approach also offers a more bio- 
logical interpretation of evolutionary dynamics than com- 
mon formulations like the Fisher- Wright or Moran pro- 
cess [TJ [31 [TH [TS] . That is to say, fitter individuals prevail 
due to higher birth rates and not by winning a tooth-and- 
claw struggle where the birth of one individual directly 
results in the death of another one. The advantage of our 
formulation is illustrated by the dilemma of cooperation 
where a transient increase in cooperation can be found 
(which does not exist in standard approaches, Eq. (lj). 

In the following, we consider two different traits, A and 
B, in a well-mixed population, all the same generalizing 
the model to more traits is straightforward. The state of 
the population is then described by the total number of 
individuals N = Na + Nb and the fraction of one trait 
within the population, x — Na/N. The stochastic evo- 
lutionary dynamics is fully specified by stochastic birth 
and death events with rates 



T^ s = Gs(x,N)Ns 



= D s (x,N)N s , (3) 



where Gs(x,N) and D$(x,N) are per capita reproduc- 
tion and death rates for an individual of type S S {^4, B}, 
respectively. We consider these rates to be separable into 
a global and relative part, meaning a trait-independent 
and trait-dependent part: 



G s =g(x,N)f s (x), D s = d(x,N)w s (x) 



(4) 
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The global population fitness, g{x, N), and the global pop- 
ulation weakness, d(x,N), affect the population dynam- 
ics of all traits in the same manner. For example, they 
account for constraints imposed by limited resources or 
how one strategy impacts the whole population. In con- 
trast, the relative fitness, fs(x), and the relative weak- 
ness, ws(x), characterize the relative advantage of one 
strategy compared to the other. They are different for 
each trait and depend, in a first approach, only on the 
relative abundance x [23) . The relative fitness terms, 
fs(x), affect the corresponding birth rates, and the rela- 
tive weakness functions, ws(x), describe the chances for 
survival of distinct traits. 

While in evolutionary game theory only the relative fit- 
ness is considered [5] , and common models of population 
dynamics take only the global functions into account, we 
here consider both global and relative fitness and show 
how their interplay determines the evolutionary outcome 
of a system. In the following, we set wa(x) = wb(x) = 1 
in order to compare our unifying approach with standard 
formulations [2] . Though the full stochastic dynamics are 
given by a master equation, it is instructive to disregard 
fluctuations for now and examine the corresponding set 
of deterministic rate equations: 

d t x = g(x,N)(f A (x)-(f))x, (5a) 
d t N = [g(x, N) (/) - d(x, N)] N, (5b) 

where (/) = x/a + (1 — x)/b denotes the average fit- 



ness. Eq. (5a) has the form of a replicator equation [2J. 
However, in Eq. ( 5a) there is an additional factor, namely 



the global population fitness g{x, N). This leads to a cou- 
pling of x and N whose implications we will discuss later 
on. Similarly, Eq. (5b) describing population growth is 
coupled to the internal evolution, Eq. ( 5a ) . Note that for 



frequency-independent global functions, g(x,N) = g(N) 
and d(x,N) = d(N), Eqs. ([5} resemble Eqs. and 
Only then, the deterministic dynamics reduces to the 
common scenario |12[ I13[ 115) . where a changing popu- 
lation size is immaterial to the evolutionary outcome of 
the dynamics [3]. For the full stochastic dynamics the 
strength of fluctuations scales as Vl/iV El E] and 
thereby is strongly affected by population growth. 

In more realistic settings, the global fitness and weak- 
ness functions, g(x,N) and d(x,N), can also depend on 
the relative abundance, x. This implies an interdepen- 
dence of population growth and internal evolution. In 
the following, we focus on one particular but very im- 
portant example: the dilemma of cooperation in a grow- 
ing population. There is an ongoing debate in socio- 
biology regarding how cooperation within a population 
emerges in the first place and how it is maintained in 
the long run [16] . Microbial biofilms serve as versatile 
model systems [Sl fTTHTO] . There, cooperators are produc- 
ers of a common good, usually a metabolically expensive 
biochemical product. For example, for the proteobacte- 
ria Pseudomonas aeruginosa, cooperators produce iron- 
scavenging molecules (siderophores) . Released into the 
environment these molecules strongly support the iron 



uptake of each individual in the population [15]. Coop- 
erators thereby clearly increase the global fitness of the 
population as a whole, leading to a faster growth rate 
and a higher maximum population size |19j . In such a 
setting, however, non-producers ("cheaters") have a rel- 
ative advantage over cooperators as they save the cost 
of providing the common good, e.g. the production of 
siderophores. Hence, their relative fraction is expected to 
increase within the population implying that the global 
fitness of the population declines. Surprisingly, as we 
show in the following, a coupling between growth and in- 
ternal evolution can overcome this dilemma transiently 
and the average level of cooperators can increase despite 
a disadvantage in relative fitness. 

We model the internal evolutionary dynamics by the 
prisoner's dilemma game [SJ [TBJ. Within this stan- 
dard approach, individuals are either cooperators (A) 
or cheaters (B). While cooperators provide a benefit 
b to all players at the expense of a (metabolic) cost 
c < b, a cheaters save the cost by not providing the 
benefit. The relative fitness of these traits is given by 
/a(x) = 1 + s [(b — c)x — c(l — x)) and fs(x) = 1 + sbx, 
respectively, where the frequency-independent and de- 
pendent parts are weighted by the strength of selection 
s [12] . Analyzing the prisoner's dilemma per se, defectors 
are always better off than cooperators because of their 
advantage in relative fitness, f A (x) < /s(x) [IB] . In the 
following, we choose for specificity 6 = 3 and c = 1, how- 
ever, our conclusions are independent of the exact values. 

Importantly, cooperation positively affects the whole 
population by increasing its global fitness, e.g. by pro- 
duction of a common good like siderophores. Here, we 
consider bounded population growth with a growth rate 
increasing with the coopcrator fraction x. In detail, we 
choose a x-dependent global fitness, g(x) = 1+px, and a 
iV-dependent global weakness, d(x, N) = N/K account- 
ing for limited resources. For p — 0, one obtains the well- 
known dynamics of logistic growth |20| with a carrying 
capacity K. For p > 0, the carrying capacity, K(l +px), 
depends on the fraction of cooperators. For instance, 
for P. aeruginosa [15], the iron uptake, and hence the 
birth rates, increase with a higher siderophore density 
and therefore with a higher fraction of cooperators. 

To analyze the evolutionary behavior of our model we 
performed extensive simulations of the stochastic dynam- 
ics given by the master equation determined by the birth 
and death rates, Eqs. All ensemble averages were 
performed over a set of 10 4 realizations. In Fig. [l] the 
average population size, N, and the average fraction of 
cooperators, x, are shown for different initial popula- 
tion sizes, No. The influence of a frequency-dependent 
growth on the population is twofold. First, starting in the 
regime of exponential growth, the frequency-dependent 
global fitness may cause an overshoot in the population 
size [Fig. TJa)]. Second, and more strikingly, the selec- 
tion disadvantage of cooperators can be overcome and 
a transient increase of cooperation emerges, [Fig. [ljb)]. 
It is maintained until a time t r , which we term as the 
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FIG. 1: The dilemma of cooperation in growing populations, 
(a) Average population size over time. Due to a cooperation- 
mediated growth advantage, it can show an overshoot. The 
red line corresponds to simulation results while the black line 
is obtained by evaluating Eqs. (Jfjjl. (b) The average level of 
cooperation increases transiently for times t < t c , especially 
if the initial population size is small meaning fluctuations are 
large. The parameters are given by xo = 0.5, b = 3, c = 1, 
s = 0.05, K = 100 and p = 10. iV is 4 (red line), 2 (blue line), 
and 12 (green line), respectively. The black line is obtained 
by evaluating Eqs. Q for No = 4. Cooperation times t c are 
denoted by thin lines of the corresponding color. 



cooperation time. 

Both phenomena rely on a subtle interplay between in- 
ternal evolution, with a selection pressure towards more 
defectors, and population growth, with a growth rate in- 
creasing with the fraction of cooperators. While the over- 
shoot in population size can already be understood on the 
basis of the rate equations, 

dtx = — s(l + px)x(l — x), (6a) 
d t N = [(1 + px) (f) - N/K] N, (6b) 

the transient increase of cooperation is a genuinely 
stochastic event as discussed in detail below. A first 
impression of the antagonism between selection pressure 
and growth can already be obtained by examining the 
characteristic time scales. While the fraction of coopera- 
tors changes on a time scale t x oc 1/s, the population size 
evolves on a time scale tm oc 1. Hence, the strength of se- 
lection, s, regulates the competition between population 
growth and internal dynamics. For s » 1, selection is 
much faster than growth dynamics. Therefore, the rapid 
ensuing extinction of cooperators cannot be compensated 
for by the growth advantage of populations with a larger 
fraction of cooperators. In contrast, in the limit of weak 
selection (s <C 1), growth dynamics dominates selection 
and both an overshoot in the population size and a tran- 
sient increase of cooperation become possible (see below) . 
In the following we focus on this latter, more interesting, 
scenario of weak selection (tjv < t x ). 

Let us first consider the overshoot in the population 
size [Fig. [TJa)]. It is caused by a growth rate and a car- 
rying capacity which are increasing functions of the frac- 
tion of cooperators (here we use p = 10 as observed in mi- 
crobial experiments |18jV For t < t x , a small population 



(N <C K(l+pxo)) with an initial fraction of cooperators, 
Xo, grows exponentially towards its comparatively large 
carrying capacity K(l+pxo). During this initial time pe- 
riod the fraction of cooperators evolves only slowly and 
can be considered as constant. On a longer time scale, 
t > t x , however, selection pressure drives the fraction of 
cooperators substantially below its initial value Xq, lead- 
ing to a smaller carrying capacity, K(\ + px). Finally, 
cooperators go extinct and the population size decreases 
to K . This functional form of N(t) is well described by 
the rate equations see black line in Fig. [l|a) . 

In contrast, the transient increase of cooperation, cf. 
Fig. ljb), cannot be understood on the basis of a sim- 
ple deterministic approach, where dtx < holds strictly 
(see black line in Fig. |TJb)). It is a genuinely stochas- 
tic effect, which relies on the amplification of stochas- 
tic fluctuations generated during the initial phase of the 
dynamics where the population is still small. In more 
detail, for small populations , the fraction of coopera- 
tors is subject to strong fluctuations and differs signif- 
icantly from one realization to another. Crucially, due 
to the coupling between the growth of a population and 
its internal composition, these fluctuations are amplified 
asymmetrically favoring a more cooperative population, 
i.e. growth, set by the global fitness g(x), is amplified 
by an additional cooperator while it is hampered by an 
additional defector. This implies that the ensemble of re- 
alizations becomes strongly skewed towards realizations 
with more cooperators. If this effect is strong enough the 
ensemble average x{t) = J2i ^A,i(t)/J2i ^i(*)> which de- 
scribes the mean fraction of cooperators when averaging 
over different realizations i, increases with time. Due 
to a subsequent antagonism between selection pressure 
towards more defectors and asymmetric exponential am- 
plification of fluctuations during growth phase, there is 
only a transient increase of cooperation in a finite time 
window, t G [0,t c ]. These findings are illustrated in a 
movie |24j showing the time evolution of the probability 
distribution for an ensemble of stochastic realizations. 

Additional qualitative and quantitative insights can be 
gained from analytic calculations via a van Kampen ap- 
proximation PP, see EPAPS document [25 . Thereby 
starting with a master equation given by Eqs. (|3) first 
and higher moments of the fluctuations can be obtained. 
They show that fluctuations during the first generation 
(i.e. doubling the initial population size on average) are 
by far the dominant source for the variance in the com- 
position of the population. In addition (see below), these 
calculations give a strictly lower bound on the parame- 
ter regime where the cooperation time is finite and thus 
quantify the magnitude of fluctuations necessary to over- 
come the strength of selection acting against cooperators. 

Fig. [2]shows the cooperation time, t c , with varying se- 
lection strength, s, and initial population size, Nq. For 
large s and No (light grey area), t c is identical to zero, i.e. 
the fraction of cooperators always decreases as predicted 
by the deterministic replicator dynamics, Eq. i(5a; In con- 
trast, if s and Nq are sufficiently small, t c is finite. The 
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FIG. 2: Dependence of the cooperation time t c on the strength 
of selection s and the initial population size No- There exist 
two distinct phases: the phase of transient maintained cooper- 
ation (where t c > holds) and the phase of extinction of coop- 
eration (where t c — 0). The boundary of both phases (solid 
line) is approximately given by sNo w p/(l + P%o) (dashed 
line). The cooperation time i c is shown for varying s but 
fixed No in the inset. See text and [25] 

transition between these regimes is discontinuous marked 
by a steep drop in the cooperation time from a finite value 
to zero; see Fig. [2] (inset). A strictly lower bound for the 
phase boundary (Fig. [2] solid line) can be derived an- 
alytically by comparing the antagonistic effects of drift 
and fluctuations, see [25]. Its asymptotic behavior for 
large N is given by sN p/(l + pxo) (Fig. [2] dashed 
line). This behavior resembles the condition for neutral 
evolution E]- Indeed, for sN < p/(l + px n ), fluctu- 
ations dominate and the system evolves neutrally. It is 
this neutral evolution leading to sufficiently large fluc- 
tuations which in turn - by asymmetric amplification - 



result in a transient increase of cooperation. 

In summary, we introduced a general approach, which 
couples the internal evolution of a population to its 
growth dynamics. Both processes originate from birth 
and death events and are therefore naturally described 
by a unifying stochastic model. The standard formu- 
lations of evolutionary game theory and population dy- 
namics emerge as special cases. Importantly, by includ- 
ing the coupling, our model offers the opportunity to in- 
vestigate a broad range of phenomena which cannot be 
studied by standard approaches. We have demonstrated 
this for the prisoners dilemma in growing populations. 
Here, a transient regime of increasing cooperation can 
emerge by a fluctuation- induced effect. For this effect, 
the positive correlation between global population fitness 
and the level of cooperation is essential. Similar to the 
Luria-Delbriick experiment |22j , initial fluctuations in the 
fraction of cooperators are exponentially amplified. Here, 
this renders it possible for cooperators to overcome the 
selection advantage of defectors. 

In biological settings, growth is ubiquitous: popula- 
tions regularly explore new habitats, or almost go ex- 
tinct by external catastrophes and rebuild afterwards. 
For a realistic description, it is therefore necessary to 
relax the assumption of a decoupled population size. Es- 
pecially for bacterial populations undergoing a life-cycle 
with a repeated change between dispersal and maturation 
phases [H IT7HT9] , a transient increase in cooperation may 
be sufficient to overcome the dilemma of cooperation. 
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edged. 
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Evolutionary game theory in growing populations 

Anna Melbinger, Jonas Cremer, and Erwin Frey 
Supplementary EPAPS document: conditions for the transient increase of cooperation 

The transient increase of cooperation emerges if initial fluctuations in the evolutionary dynamics are sufficiently 
large such that the asymmetrical amplification of those can overcome the selection advantage of cheaters. In this 
Supplementary Material we derive the conditions for the transient increase. In particular, we give an analytical 
expression for the phase boundary in Fig. 2 (black line). 

The full stochastic dynamics is given by the master equation determined by the birth and death rates, Eq. (3), 

T $ ^ A (A-l,B)(A-l)P(A-l,B) + r $ ^ B (A,B-l)(B-l)P(A,B-l) 

+ T A ^ $ (A+1, B)(A + 1)P(A+1, B) + T B ^(A, B+1)(B+1)P{A, B + l) 
- [T^ A (A, B)A + T^ B (A, B)B + T A ^(A, B)A + T B ^(A, B)B] P(A, B). 

(7) 

Here, A = N A and B = N B stand for the number of individuals of both traits. We approximate the master equation 
upon performing a van Kampen expansion pQ. To this end, we consider A and B as extensive variables which we 
write as 

A = n a (t) + Vn£ , 

B = Clb(t) + Vttfi . (8) 



dP(A,B) 
dt 



Here, f2 is of the order of the actual system size, and deterministically evolving densities a(t) and b(t) are corrected 
by fluctuations £ and /i. By this Ansatz the strength of fluctuations is correctly considered; their relative impact 
decreases like with increasing system size. In the following, we consider the initial dynamics of the population 

when starting with a small population size No. Then, £1 is of the order £1 ss Nq. Death events can be neglected as 
the initial population size is far below the carrying capacity, Nq/K w 0. 

To proceed, we expand Eq. (7) in orders of The deterministic equations follow to leading order, O 

see Eqs. (6) with N/K -> and x(t) = a(t)/[a(t) + b(t)}. The next leading order, 0(f2°), results in a Fokker-Planck 
equation for the probability distribution of the fluctuations, n(£,/i). The dynamics in n(£,/i) is coupled to the 

deterministic equations and can be extended to include higher orders, O (l/VTlj . From the Fokker-Planck equation 

for n(£, fj,), differential equations for the first moments of £ and fi can be obtained. They have the following functional 
form, 



$1 



d t (fi) =0,(0 + D 2 (fx) + -^=(D 3 (e) + £>4<£M) + £ 5 (M 2 » + O 



(9) 



The constants Cj and Di with i G {1,2,3,4,5}, depend on the parameters s , 6, c, p, the deterministic parts of the 
composition of the population, x(t) = a(t)/[a(t) +b(t)], and the population size n(t) = a(t) + b(t) (in units of O), 

respectively. Importantly, the second moments couple into the dynamics only through O corrections. 

Neglecting these second and higher order moments, the ensuing linear equation has an unstable fixed point at 
((£), (//))* = (0, 0). The eigendirection with the larger (positive) eigenvalue has a component in the ^-direction which 
is significantly larger than its component in the /x-direction. As a consequence, the fluctuations in the number of 
cooperators (£) are amplified more strongly than those of the defectors (//); fluctuations are asymmetrically amplified. 

Next, we consider the effect of the second moments on the dynamics. Consider a single initial state without any 
variance (and all other higher moments identically zero), starting the dynamics in the fixed point, ((£), (fJ>))* = (0, 0). 
Then, since the first moments are zero, only higher orders in Eq. ^ lead to deviations from the (linearly unstable) 
fixed point. Once such deviations are generated these are amplified by the (linearly) unstable dynamics, i.e. the first 
moments in Eq. ([9|. In more detail, consider the differential equations of the second moments which, for t — > 0, have 
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the following asymptotic form: 

<9 t (£ 2 ) =n(l + px) [l + s(bx-c)]x, 

d t m =0, 

d t (n 2 ) =n(l+px){l + sbx)(l - x). (10) 

Starting with zero at t — 0, both, (£ 2 ) and (li 2 ) increase linearly in time (note that the fitness of a cooperator l + s(6:r — 
c) > since otherwise the birth rate would be negative). Within one generation, t g = 1/ [(1 + px)(l + s(b — c)x)] 
(compare Eq. (6b)), i.e. doubling the population size on average, finite variances (£ 2 ) ff and (Li 2 } g are generated. This 
variance can be taken as a lower bound. We even expect it to be a reasonable estimate for the actual value since the 
impact of the variance created in following generations on Eqs. [9] is strongly suppressed by the increase in population 
size. 

Upon inserting the values (£ 2 ) g and (n 2 ) g into Eq. ^ one can now calculate the time evolution of the first moments, 
(£) and (li). This allows to determine the conditions necessary for a transient increase of cooperation by analyzing 
the fraction of cooperators ( A + B ); see Eqs. (j8j). The phase boundary separating the regimes of transient increase and 
immediate decrease of cooperation is defined by an initially stationary fraction of cooperators: dt ( A + B ) = at t w 0. 

The ensuing phase boundary is plotted in Fig. 2 (black line). The deviation from the actual (numerically determined) 
transition line is small for intermediate Q and goes to zero for larger Q. By evaluating the expression in orders of s/p, 
the lower bound of the transition line can be further simplified. To first order one finds 

° = m^y <"> 

with fin = No; see Fig. 2, dashed line. Note that this expression gives the asymptotically correct results for large f2. 

It is instructive to compare this result with the theory of neutral evolution [5] where a condition sN oc 1 separates 
regimes of neutral and selection-dominated evolution In the present case, for the transient increase of cooperation 

to occur, the system has to evolve neutrally in the initial phase to create a large enough variation in the fraction 
of cooperators. Then, after being asymmetriclly amplified, these fluctuations can overcome the selection pressure 
towards more defectors. This is mathematically reflected in Eqs. ([9| and (10 1. Initially, the second moments increase, 
Eqs. (10), which then feeds into Eqs. Q and lead to an increase in the first moments. Finally, the good agreement 
of the phase boundary with its lower bound, reassures that the variation in cooperators fraction is mainly generated 
at the beginning of the dynamics. 
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